Fast density-matrix based partitioning of the energy over the atoms in a 
molecule consistent with the Hirshfeld-I partitioning of the electron density 

Diederik Vanfleteren/'^'[^ Dieter Ghillemijn,'-^ Dimitri Van Neck/'^ 
Patrick Bultinck/'^ Michel Waroquier/'^ and Paul W. Ayers"* 

^Members of the Ghent Brussels Quantum Chemistry and Molecular Modeling alliance. 
^ Ghent University, Center for Molecular Modeling, 
Technologiepark 903, B-9052 Zwijnaarde, Belgium 
^ Ghent University, Department of Inorganic and Physical Chemistry, 
Krijgslaan 281 (S3), B-9000 Gent, Belgium 
'^McMaster University, Department of Chemistry, 
Hamilton, Ontario L8S 4M1, Canada 

Abstract 

For the Hirshfeld-I atom-in-molecule model, associated single-atom energies and interaction energies at 
the Hartree-Fock level are determined efficiently in one-electron Hilbert space. In contrast to most other 
approaches, the energy terms are fully consistent with the partitioning of the underlying one-electron den- 
sity matrix. Starting from the Hirshfeld-I atom-in-molecule model for the electron density, the molecular 
one-electron density matrix is partitioned with a previously introduced double-atom scheme [Vanfleteren D. 
et al., J Chem Phys 2010, 132, 1641 1 1]. Single-atom density matrices are constructed from the atomic and 
bond contributions of the double-atom scheme. Since the Hartree-Fock energy can be expressed solely in 
terms of the one-electron density matrix, the partitioning of the latter over the atoms in the molecule leads 
naturally to a corresponding partitioning of the Hartree-Fock energy. When the size of the molecule or the 
molecular basis set does not grow too large, the method shows considerable computational advantages com- 
pared to other approaches that require cumbersome numerical integration of the molecular energy integrals 
weighted by atomic weight functions. 
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I. INTRODUCTION 



In the last decades, several algorithms have been developed to identify the atom in the molecule 
(AIM). For example the MuUiken method ([U relies on the attachment of basis functions to atomic 
centers; Natural Population Analysis [2] relies on the analysis of blocks of the one-electron density 
matrix (IDM) expressed in some molecular orbital basis; some other methods rely on the parti- 
tioning of the molecular electron density in AIM parts. However, not all AIM properties can be 
directly expressed in terms of the electron density. Of course the density determines all the prop- 
erties and there are even explicit formulas for them, but these formulas are often computationally 
impractical. A common example is the kinetic energy of an AIM, which is directly computable 
from the full IDM p (r, r'), but not (trivially) from the electron density, the diagonal element of 
the IDM. 

A widely adopted solution, which is also the most common recipe employed within the Quan- 
tum Theory of Atoms in Molecules (QTAIM) [|3jl, is to partition such molecular properties us- 
ing the same atomic weight functions wa (r) that are used to partition the molecular density. 
QTAIM uses a zero flux condition of the electron density to define interatomic surfaces, result- 
ing in nonoverlapping atomic regions. Atomic energies are obtained by partitioning the kinetic 
energy over the atomic domains and using the virial ratio to construct the corresponding atomic 
potential energy. Other approaches also rely on atomic regions to partition the molecular energy, 
but use them in a more general way to decompose the energy into one- and two-atom terms [|4]- 
13. Recently, Mandado et al. [8] partitioned the Hartree-Fock energy in terms of the overlapping 
Hirshfeld atoms. Their scheme appears useful to investigate proton acidity, the anomeric effect 
and group transferability, it has also been used in studies of bonding and polarizability |l9l [TOl . 
However, partitioning of the molecular energy in atomic fragments can be ambiguous as it often 
relies on the introduction of an arbitrary number of partitionings of unity [jTTl| into the energy ex- 
pressions. Moreover, the exact place where the partition of unity is introduced in an expectation 
value expression can have an important influence on the resulting AIM condensed values ffT2l . 
This ambiguity is often circumvented by the convention to introduce the unity and its partitioning 
into weight functions (1 = YIia ^a(?')) before any operator [fTT]|. 

To avoid these problems, we partition the Hartree-Fock energy starting from a partitioning of 
Hartree-Fock molecular IDM. At the Hartree-Fock level of theory the energy can be expressed 
solely and directly in terms of a (partitioned) molecular IDM. A partitioning of the Hartree- 
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the molecular fragments", J dv KpA (r, r') 
/ dr wa (r) Ap (r, r') 



Fock molecular IDM therefore directly leads to a natural partitioning of the energy. In a pre- 
vious paper [fT3l we introduced a double-atom partitioning scheme for the molecular IDM that 
is consistent with existing partitioning schemes for the molecular density. In the current work 
we use molecular IDM fragments from that scheme to calculate the energy terms naturally as- 
sociated with these fragments. The strategy of this work is therefore to calculate "properties of 

, instead of "fragmenting molecular properties", 

r' 

This is in line with the work of Bultinck et al.[fT2| for the derivative 
of the density function with respect to a change in the number of electrons. We also examine the 
correspondences and differences between both approaches. 

From a practical point of view, a fragmentation of molecular properties requires that the ex- 
pectation value integrals are computed numerically on a large spatial grid. In general there is no 
simple analytical expression for the atomic weight functions and these weight functions are often 
not well expressed in one-electron Hilbert space. In contrast, molecular IDM fragments show a 
satisfying basis set convergence in one-electron Hilbert space ffTSl. Therefore, it is tempting to 
calculate the atomic and interaction energies of the IDM fragments in one-electron Hilbert space, 
avoiding cumbersome numerical integrations in r-space. As a consequence, we expect to find that 
significant computational advantages are a major asset of our current approach. 

n. ATOMIC DENSITY MATRICES THAT ARE CONSISTENT WITH THE HIRSHFELD-I 
PARTITIONING OF THE ELECTRON DENSITY 

We use the notation x = ra io specify the single-electron states in coordinate space, where a 
represents the spin degrees of freedom. The one-electron density matrix (IDM) for an A^-electron 
molecule with wave function ^!ixi, . . . .x^) h defined as 



We restrict ourselves to molecules with a singlet ground state. In that case 



(1) 



p{x,x) = ]^5^^„^p{r,r'), (2) 

and the electron spin can be discarded. In a previous paper we introduced a double-atom partition- 
ing scheme for the molecular spin- summed IDM [fT3ll . 

p{r,r') = ^PABir,r'), (3) 

AB 



in terais of atomic {A = B) and diatomic contributions {A ^ B), where A and B label atoms. The 
individual contributions are defined as 

pAB{r, r') = ^ [wA{r)wBir') + WB{r)wA{r')] p(r, r'), (4) 

with positive local weight functions wa(^) obeying 

J2^A{r) = l. (5) 

A 

At the Hartree-Fock level the energy is expressed solely in terms of the IDM. Therefore, the 
partitioning of the latter over the atoms in the molecule leads naturally to a corresponding parti- 
tioning of the Hartree-Fock energy. However, for a simple energy partitioning we require that (1) 
it is based on single-atom density matrices (rather than double-atom density matrices) and (2) the 
electron density of the atoms is local and positive definite. Note that any single-atom density ma- 
trix pAir, r') unavoidably has bad localization properties (in contrast to the double-atom density 
matrices) [[T3l . and therefore we restrict the requirement of locality and positive definiteness to 
the electron density, the diagonal element of the IDM. The latter requirement is needed to prevent 
the energy components to become unrealistically large on the chemical energy scale. When the 
single-atom density matrices are defined as 

PA{r,r') = '^pAB{r,r'), (6) 

B 

it is clear that on the diagonal (r = r') the following definition is obtained, 

PAir) = pAB{r, r) = WA{r)p{r) (7) 

B 

that is familiar from existing partitioning schemes for the electron density (with good localization 
properties). Also the fundamental property 

Y,PA{r,r') = p{r,r') (8) 

A 

is obeyed as a trivial consequence of Eq. ([3]). Note that these single atom density matrices may 
not be N-representable, an issue that, however, does not stand in the way of obtaining integrated 
quantities such as energy contributions. 

In order to define a computationally efficient energy partitioning scheme, we assume real wave- 
functions and express the partitioned density matrices in the finite basis set used for the molecular 
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calculation. E.g., in the molecular Hartree-Fock basis set, the following spatial integrals 

Sf^ = j dr<P,ir)wAir)<Pjir), (9) 

that represent elements of the regular atomic overlap matrix (AOM), are sufficient to determine 
the single-atom IDM in this basis, 

iPAh = j drdr'Ur)pA{r,r')<P,{r'). (10) 

The precise expressions for the density matrices are: 

{Pa\, = J2^-lj2^^^nSfn + Sf^Sfj = ^-{d, + d,)Sf^, (11) 
n B 

where a simple sum rule 



is used to simplify the expression in Eq. (11). 



In order to ensure that the single-atom density pAir) is local, it can be made to coincide with 
the AIM density of some well-established density partitioning scheme like Hirshfeld [fT4l . iterative 
Hirshfeld [Ull, Iterated Stockholder Atoms IMO or QTAIM [IlIIHlIia, 

^(r) = hA{r)pir) (13) 

in which the AIM density is obtained by multiplying the molecular density p(r) with the charac- 
teristic weight function for that AIM technique hA{r). One simply takes 

WA{.r) = hA{r). (14) 

An alternative to Eq. Q, consisting of distributing the diatomic contributions in a weighted man- 
ner, was also investigated in [13], but will be discarded here because it gave inferior results. 

III. ENERGY DECOMPOSITION 

A. Self-energies and interaction energies of the IDM fragments 

Following the ideas behind the Interacting Quantum Atoms (IQA) [|20] - l22l . but from the per- 
spective of density matrix fragments in one-electron Hilbert space, each single-atom IDM with 
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elements {pA)ij constructed in section |n] can be thought of as forming an atomic subsystem [|23ll 
when combined with the nucleus on center A. The molecular energy can then be decomposed as 
the sum of the "self-energies" of the atomic subsystems and the interaction energies between them: 

^moi ^ J2 Ef^ + ^ E%. (15) 

A A<B 

At the closed-shell Hartree-Fock level, the self-energy of such an atomic subsystem corresponds 
to: 

Eff = ^(t + VX^'UpAh + ^ E ^^1^^ (iPAy^kipAh - lipA)uipAhk) (16) 
ij ijkl ^ 

where 

tij = j (f)i{r)i(j)j{r)dr 

J \r — ±tA\ 

Vijki = I 4>,{ri)(f)Ar2)-. ^ 74>k{ri)Mr2)dridr2 

(17) 

In Eq. (17), £ is the kinetic energy operator and Za is the nuclear charge on an atom A with 



nuclear coordinate R^- Vijki are two-electron integrals that are not anti-symmetrized. Note that it 
is implied in Eq. ( [76] ) that the atomic subsystems are spin-averaged, consistent with Eq.([2]j4]). The 
atomic subsystems interact with each other, according to the following expressions: 

' AB 

+ '^Vi.jkl UpA)ik{pB)jl - ^{pA)il{pB)jkj (18) 
ijkl ^ ^ 

^int ^ J^Eab- (19) 

A<B 

j^mt represents the total interaction energy of all subsystems within the molecule. quantifies 
all interactions between atom pairs in molecules, including the interaction between atom pairs that 
do not share a chemical bond. This quantity is also useful to assess the strength of the interactions 
in a ring structure. E'^j^ does not depend on the choice of reference atoms, it only depends on the 
AIM. The localized character of the Hirshfeld-I densities ensures that the promotion and interac- 
tion energies are within a reasonable range of values, although they are often significantly larger 
(in absolute value) than typical "bond energies". 



B. Promotion energies 



Since it is well known that the molecular environment induces only relatively small changes 
in atomic energy [j23|, atomic energies are usually referenced to the energy of the isolated atoms. 
Atoms within a molecular environment are slightly distorted compared to the isolated atoms. This 
implies that their energy, with respect to the Hamiltonian of the isolated atoms, has increased. 
This energy increase is called the atomic promotion energy E^°"^. When the atom in a molecule is 
identified with the atomic subsystem defined above, its promotion energy is the difference between 



the (Hartree-Fock) atomic self-energy E^''^ in Eq. (16) and the Hartree-Fock energy of the 
isolated atom: 

^prom ^ ^self _ ^0 (20) 

^prom ^ (21) 

A 

The atomic promotion energies can be considered to result from three successive processes: 
first a charge transfer step that accounts for the fact that the AIM has a different atomic charge 
from the isolated atom. 

AE^^ = EliQj,) - E\{Q) ill) 

Classically, the reference for this step is the neutral isolated atom in its ground state -E'°(0). 
E\{Qa) represents the charged isolated atom, and is approximated as the linearly interpolated 
value between the energies of the isolated atoms with an integer number of electrons N < {Qa = 
N + a) < + 1. In principle, the HF energy as a function of the number of electrons N is a 
concave curve between the integers [I8l l24l - l27ll . The assumption that it is linear is based on the fact 
that this holds for the exact energies [|28l - [30ll . This leads to E^{Qa) computed as: 

E'HQa) = aEl{N + 1) + (1 - a)El{N). (23) 

The interpolated state E\{Qa) is characterized by a value of (5*^) (where S is the spin angular 
momentum) that does not always correspond to a singlet whereas the AIM is always assumed to 
be spin averaged (see Eq. Q). We therefore introduce a second step that accounts for the spin 
averaging. It corresponds to the energy difference between on the one hand the isolated atom with 
interpolated charge and averaged spin E^jf{QA) and on the other hand E\{Qa), 

AeJ = EfiQA) - E\{Qa). (24) 
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E'^^{Qa) is obtained as the interpolated value between the energies of the spin-averaged isolated 
atoms with an integer number of electrons, 

EfiQA) = aEfiN + 1) + (1 - a)EY{N), (25) 
and the E'^^{N) are calculated from the following formula: 

A ij ijkl ^ ' 

where pA is the spin-summed IDM of the isolated atom with N electrons. 

In the third step, charge reorganization takes place which corresponds to the difference between 
the self- energy of the of the atomic IDM fragment and the energy of the isolated atom with 
interpolated charge and averaged spin, 

AE^^ = Eff-Ef{QA). (27) 
Finally, the promotion energy is retrieved as the sum of three terms as shown below: 

^pjom ^ ^^CT ^ ^^J ^ ^^CR_ (28) 

Although one may point out that there is some arbitrariness in the order chosen for the processes. 



the decomposition in Eq. (28 1 might give a general idea about their relative importance. 



C. Bond energies 

Since the promotion and interaction energies depend on the AIM definition, they are not directly 
observable. The sum of all promotion and interaction energies, i.e. the atomization energy AE'^^ 
of the molecule, can be compared to experiment: 

^^at ^ ^moi -^El = Y^ E^"^ + J2 Eab- (29) 

A A A<B 

In order to get a quantity that is more in line with the chemical concept of a "bond energy", it 
would be convenient to interpret the atomization of a molecule strictly in terms of atom pairs. 
Therefore, the atomic promotion energies {E^°"^) are attributed to the atom pairs (A < B), 

/ rpint \ / rpint \ 

j^^^ ^ j^,rom ^ j ^ j^,rom (^^-^ ] , (30) 



and Hartree-Fock bond energies are derived, 

^AB = ^AB + ^AB ) (31) 

that reproduce exactly the Hartree-Fock atomization energy 

Ai?("^^) = E^". (32) 

A<B 

These bond energies can be compared to average bond dissociation energies. 

rv. CONSISTENT DECOMPOSITION OF MOLECULAR PROPERTIES 

In the introduction we already noted that the current approach to calculate self-energies and in- 
teraction energies for IDM fragments is unambiguous. That is not the case for the widely adopted 
approach to partition the molecular energy using numerical integration of the molecular integrals 
weighted by the atomic weight functions. In this section we explore the exemplary cases of two 



components of the atomic self-energy (see Eq. (16)): the atom-condensed kinetic energy and the 
atom-condensed Fock energy. 

A. Atom-condensed kinetic energy 

The strategy adopted in the current paper is to calculate energy components of the molecular 
fragments, rather than to fragment the molecular energy. For the kinetic energy in particular, this 
can have important consequences. The kinetic energy t^ of the single-index atomic density matrix 
is calculated in Eqs. ([T6f [17]) in a finite basis set, but the corresponding expression in r-space is: 

^A = -\j cirVV(^,OL=,,, (33) 

where the notation 1^=^/ indicates that r' is replaced by r after the action of V^^-, on Pa(^, f') 
but before the integration is carried out. Note that different but equivalent representations of the 
kinetic energy operator exist, i.e. 

tA=\j t/r V-VW(r,r')|,=,,, (34) 



should give the same result as in Eq. (33). This is indeed fulfilled for the present formulation 



in terms of atomic density matrices, as follows directly from partial integration by generalizing a 



(35) 



well known relationship for the kinetic energy of the molecular IDM, 

to its fragments. In contrast, the expression for the fragmentation of the molecular kinetic en- 
ergy with Hirshfeld-I weight functions clearly depends on the representation of the kinetic energy 
operator [[3111321, since in general 



yields a result that differs from 

ti = \j dr WA{r)V-Vp{ry)l^^,. 



(36) 



(37) 



Only in special cases, e.g. when QTAIM weight functions are used, do Eq. ( [36| ) and Eq. ( |37| ) 
coincide. This shows that a naive fragmentation of the molecular energy by introducing weight 
functions in the molecular integrals is ambiguous by nature, in contrast to the fragmentation of the 
molecular IDM and the calculation of the associated energy components. 



It is shown in the appendix that the r-space expressions for Ia (Eqs. (33) and (34) ) and 



(Eq. (36) ) are mathematically identical. Since Eq. (36) does not contain the action of the kinetic 
energy operator on an atomic weight function it can be calculated quite easily using 3D numerical 
integration on a spatial grid. This provides an indication of the error induced by using the finite 



basis set expression in Eq. ([T6|) instead of the full r-space expression in Eq. (33 ). 



B. Atom-condensed Fock energy 

The Fock operator provides another example of a nonlocal operator that causes problems in the 
common approach, where atomic weight functions are inserted in the molecular Fock integrals. 
Indeed, there are several possibilities for the fragmentation of the Fock energy, 

j^Fock = _1 /■ drdr' MZIlI^. (38) 
2 J \r-r'\ ^ ' 

Depending on the place where the decomposition of unity is inserted, one could have 



^Fock = J2-l[ drdr' WA{r)wB{r') 
AB ^ ^ 

= y.-] I drdr' [wA{r)wB{r) + WA{r')wB{r')] ""J' 
^-^ 4 / \r — r' 

A D II 



(39) 
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Only one fragmentation is consistent with the underlying partitioning in Eq. ([8]) of the molecular 
IDM into single-index atomic density matrices, 

1 f pA{r,r')pB{r,r') 



AB ^ 



y^.~l [ [wA{r)wB{r) +WA{r')wB{r') 



+ WA{r)wB{r') + WB{r)wA{r')] "^P^J^. (40) 

\r — r'\ 

Note that the Fock contributions with indices (A = B) are attributed to E'^''^ , while the Fock 
contributions with indices (A ^ B) are attributed to ^'ab- 



V. A FAST DECOMPOSITION IN HILBERT SPACE 

In a previous paper [fT3l . it was observed that the molecular IDM fragments show a satis- 
factory basis set convergence when they are expressed in one-electron Hilbert space. To avoid 
cumbersome numerical integrations in r-space, atomic self-energies and interaction energies are 



calculated efficiently in one-electron Hilbert space (see Eqs. (16) and (18)). The widespread al- 
ternative strategy to fragment the molecular energy (by inserting atomic weight functions in the 
molecular expressions) requires that the atom-condensed energy integrals are computed numeri- 
cally on a large spatial grid. In this section we explore the computational consequences for two 
components of the atomic self-energy (see Eq. ([16])): the atom-condensed kinetic energy and the 
atom-condensed Fock energy. 

A. Atom-condensed Kinetic energy 

When the kinetic energy is calculated by numerical integration of the molecular energy integrals 



weighted by the atomic weight functions (see Eq. 36 ), 



(41) 



it is clear that only few integrals have to be computed since nonzero contributions are obtained 
only for occupied MO's. This is an important computational advantage over using an atomic 



density matrix as in Eq. ( 16 17 ), to which all orbitals in the basis contribute. However, experience 



shows that the numerical computation of the integrals above requires quite large integration grids 
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whereas the atomic overlap integrals can be computed with sufficient precision using more modest 
grids. As will be shown, there is a trade-off between either computing numerically fewer integrals 
that require larger grids and more atomic overlap integrals that are individually computed much 
easier. 

The strategy adopted in the current paper is to calculate energy components via matrix ma- 



nipulations in one-electron Hilbert space (See Eqs. (16) and (18) ). Using Eq. (11 ), the precise 
expression for the kinetic energy is, 

= I + dj)tijSfj. (42) 

ij 

A non-zero contribution is obtained if either di or dj is nonzero. This implies that the summation 
essentially runs over all i and j and that all of the following integrals must be computed: 



Sf. = J dr<f)i{r)wA{r)(l)j{r) 
= -IJ <l>,ir) {V%ir)) dr 
However, the Uj can be calculated analytically and moderate grids suffice to construct the S:^. 



B. Atom-condensed Fock energy 

For the partitioning of e.g. the Fock-energy, the considerations of the previous subsection are 
slightly more outspoken than for the kinetic energy. Using the strategy of numerical integration, 

Eab' = / [wA{r)wBir) + WAir')wBir') 

+ WA{r)wB{r ) + WB{r)wA(r )\ 17377| — ' ^^^^ 

the number of integrals that has to be computed is limited to the square of the number of occupied 
orbitals. On the other hand, the cost of numerical integration is squared with respect to the situation 
for the kinetic energy since the integrals run over r and r' . 

When the Fock energy is partitioned via matrix manipulations in one-electron Hilbert space, 
the precise expression is: 

Elt = E(^^ + ^^)(^^- + dk)V.jkiS^Sf„ (44) 

ijkl 
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where Vijki is calculated analytically. The summation runs over all ijkl, but the 5*^^ are calculated 
on a moderate grid that runs over only r. Except for large systems, the approach presented in this 
paper is more appealing from a computational point of view. Note that E^^'^*^ is part of E'/^ when 
(A = B), while it is part of E'^^ when (A ^ B). 



VI. COMPUTATIONAL METHODS 



The formalism described in the sections |Il III was applied to the set of 25 molecules listed in 
table |l} representing a variety of chemical bonds. All molecules have a singlet ground state, apart 
from O2, for which we consider the singlet (excited) state. 

The molecular IDM was calculated at the restricted open-shell Hartree-Fock (ROHF) ll33] - [35ll 
level of theory using the (cartesian) Aug-cc-pVTZ basis, including a geometry optimization. All 
molecular one-electron density matrices were calculated with the Gaussian 03 program [|39l . 



H2 


N2 


HCl 


CH3OH 


H2O 


F2 


LiH 


CH4 


H2CO 


(H20)2 


C12 


LiF 


CH3CH3 


CHOOH 


H2O2 


Lis 


NaCl 


CH2CH2 


CO2 


NH3 


O2 


HF 


CHCH 


CO 


N2H4 



TABLE I: List of molecules in the test set. 



The partitioning scheme was implemented using the atomic weight functions hA{r) in Eq.( 13 ) 
from a Hirshfeld-I analysis. To be consistent with our previous work on the Hirshfeld-I compatible 
IDM partitioning lfT3l . the iterative Hirshfeld weights hA{r) and the Sf^ coefficients derived from 
these weights were calculated on atom-centered grids, using a logarithmic radial grid of 100 points 
with r^in = 10~^ A and r^ax = 20 A , and the 170-point Lebedev angular grids, [|40l - |45l with each 
radial shell given a randomized orientation. The sum rule of Eq.([T2]) was fulfilled with a reasonable 
precision of lO^'^ except for some very diffuse virtual orbitals (as we work in an augmented basis 
set), where the sum rule is fulfilled with less precision (10^^) because the grid is too localized. 
However, as these orbitals are of little relevance in the calculation of the energies, this does not 
significantly affect the outcome of our analysis. In order to reproduce the molecular energy in Eq. 



fT5|) exactly, the sum rule of Eq.(12) was enforced by a renormalisation of the atom-condensed 
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overlap matrices. 



VII. RESULTS 



A. Energy decomposition using the single-atom density matrices 



Table |Il] displays the total promotion and interaction energies belonging to the (single-index) 
density matrix fragments (see section III). 





r'pro'm 
'^{HF) 


^(HF) 


'-^^(HF) 






'^{HF) 


pint 
^(HF) 


^^(HF) 


-AE?^^ 
'-^^cc 


H2 


0.18 


-0.31 


-0.13 


-0.17 


CH2-CH2 


2.24 


-2.93 


-0.69 


-0.88 


F2 


0.75 


-0.71 


0.04 


-0.05 


CHCH 


1.84 


-2.32 


-0.48 


-0.63 


C12 


0.56 


-0.60 


-0.04 


-0.09 


CH3OH 


2.25 


-2.86 


-0.60 


-0.80 


Li2 


0.11 


-0.12 


-0.01 


-0.04 


H2CO 


1.85 


-2.28 


-0.42 


-0.58 


02 


1.37 


-1.35 


0.01 


-0.13 


CHOOH 


2.94 


-3.48 


-0.55 


-0.78 


N2 


1.60 


-1.79 


-0.19 


-0.34 


CO2 


2.61 


-3.02 


-0.41 


-0.60 


LiH 


0.35 


-0.40 


-0.05 


-0.09 


CO 


1.29 


-1.58 


-0.29 


-0.40 


LiF 


0.31 


-0.46 


-0.15 


-0.22 


H2O 


1.06 


-1.32 


-0.26 


-0.36 


NaCl 


0.05 


-0.16 


-0.12 


-0.15 


(H20)2 


2.16 


-2.68 


-0.52 


-0.73 


HF 


0.45 


-0.61 


-0.16 


-0.22 


H202 


1.79 


-2.02 


-0.23 


-0.41 


HCl 


0.35 


-0.48 


-0.13 


-0.17 


NH3 


1.48 


-1.80 


-0.32 


-0.46 


CH4 


1.39 


-1.91 


-0.53 


-0.66 


N2H4 


2.72 


-3.16 


-0.43 


-0.68 


CH3CH3 


2.61 


-3.50 


-0.89 


-1.12 













TABLE II: Total promotion and interaction energies (E^^^ and E*^'^^^ calculated at the ROHF/Aug-cc- 
pVTZ level of theory. For comparison with the HF atomization energies AE'^^^^ , the CCSD(T) 
atomization energies ^Ef^^ are also presented. All energies are in Hartree. 



The sum of the total promotion and interaction energies yields (minus) the Hartree-Fock at- 
omization energy -AE^j^py in most cases this represents about 3/4 of the atomization energy 
calculated at the CCSD(T) P6] - [5T1l level of theory -AE'^^. Note that in some cases the Hartree- 
Fock atomization energy is far from satisfying, e.g. it predicts that F2 has a higher energy then two 
isolated F-atoms in their ground state. The magnitude of the promotion and interaction energies is 
within the range of typical IQA (Interacting Quantum Atoms) Il20] - l22ll values, studied by Pendas 
et al. [|23l for different methods including the Hirshfeld AIM method. 
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TABLE III: Charge transfer, spin-averaging and charge redistribution energies (AE^'^, A£^^ and AE^^, 
in Hartree) for the molecular fragments of some small molecules calculated at the ROHF/Aug- 
cc-pVTZ level of theory. These are the components of the atomic promotion energies 
gpr-om rpjjg reference energy, E^, is the energy of the neutral isolated atom calculated at 
the ROHF/Aug-cc-pVTZ level of theory. Qa is the charge of the single-index atom. 
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Table |ni| displays the charge transfer (CT), spin averaging (S) and the charge redistribution 
(CR) components of the individual atomic promotion energies. The current scheme is based on a 
partitioning of the IDM. To calculate the energy components of the IDM fragments, it is necessary 
to specify the electron spin of the fragments. For singlet molecules, the current scheme averages 
the electron spin of the IDM fragments. In an attempt to avoid the spin- averaging step and get 
lower atomic promotion energies, one could implement an alternative IDM partitioning scheme 
and request that e.g. the hydrogen atoms keep their doublet structure in a H2 molecule. However, 
such requirement would lead to delocalized electron densities for the AIMs and would exclude 
consistency with the (well-localized) Hirshfeld-I model. Since the localization of the AIM densi- 
ties is a necessary condition to produce chemically meaningful results within the IQA approach 
II23II . it is clear that the spin must be averaged at least to some degree. Note that the problem re- 
lated to the spin is not addressed by the common energy decomposition schemes that are based on 
a partitioning of the molecular expectation values using AIM weight functions. In these schemes, 
the underlying electronic structure is not explicitly treated (and possibly not consistent). Energet- 
ically, it appears that the spin averaging energy defined in Eq. (|24|) is of similar importance as the 



charge redistribution energy. E.g. the average of the absolute values in Table III for AE^'^, AE'^ 



s 



and AE^^ is 0.09, 0.19 and 0.19 Hartree respectively. The charge redistribution (CR) energies 
are mostly positive, although small negative values are found for the hydrogen atoms. At first 
sight, one would expect these CR-energies to be positive for variational reasons, since it is based 
on the atomic HF-Hamiltonian, for which the isolated atom IDM with fractional charge should be 
optimal. However, is not necessarily N-representable [52], so it is variationally not sufficiently 
constrained. It is nevertheless remarkable that these negative values are rather small. 



In Table IV the interaction energies between the molecular fragments are presented. The Fock 
component is listed separately, 

FTb = E ^^^'^^ (lipAUpBhk) , (45) 

ijkl ^ ^ 

as it seems to be a robust indicator of the covalent character of that interaction. It clearly singles out 
the ionic species with a small value, whereas for the (covalent) homonuclear diatomic molecules 
in the analysis, it correlates linearly with the interaction energy (R^=0.997). In contrast to the 
SEDI (bond order) index [|53l - i62l the value is also low for the covalent but weakly bound Li2. 



The bond energies are also listed in Table IV They do depend on the choice of reference 



atoms, but chemists are more familiar with their typical magnitudes, and much chemical reasoning 
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is based on this type of quantity. The linear correlation between the interaction energies and 
the bond energies is not strong (R^=0.6). Both quantities can single out chemically not bonded 
atoms (having very small values) and for both quantities similar values are obtained for chemically 
similar bonds. The C-H bond has interaction energies between 0.45 Hartree (in CH4) and 0.40 
Hartree in (CHOOH) and bond energies between 0.12 and 0.8 Hartree. Interaction energies for 
the single, double and triple C-C bonds in ethane, ethene and ethyn are respectively -0.58, -0.98 
and -1.38 Hartree, while the corresponding bond energies are 0.10, 0.17 and 0.25 Hartree. The 
interaction energies are theoretically appealing in the sense that they follow immediately from 
the density matrix partitioning. They should not be confused with the bond energies, which are 
more in line with chemical intuition but require the introduction of isolated reference atoms. The 
analysis presented in this work is particularly useful to investigate substituent effects in molecules. 
When an atom or functional group is substituted in the molecule, the strength of all bonds in the 



molecule is affected. The bond energies defined in Eq. (31 1 are a measure for this effect. E.g. it 
is clear from the entries in Table |IV] that the C-H bond in ethene (CH2CH2) is stronger than the 
C-H bond in formaldehyde (CH2O) and formic acid (CHOOH). It is also possible to compare the 
interactions between atoms that have no chemical bond in the molecule. Although we recognize 
that it is possible to get similar information from other schemes that were developed to partition the 
Hartree-Fock energy, we stress the consistency of this approach, where different energy terms are 
derived from the same, symmetrical, atom condensed density matrices. Without this mathematical 
consistency, the results are necessarily ambiguous. 
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TABLE IV: Interaction energies (E^'^, in Hartree) for the molecular fragments of some small molecules 
calculated at the ROHF/Aug-cc-pVTZ level of theory. In the column labeled "bond", the 
notations "-", " " and "• • • " indicate respectively a chemical bond, a hydrogen bond and the 
interaction between a pair of non-bonded atoms. For comparison, the Hirshfeld-I SEDI-index 
is included as a bond order indicator and the Pock part of the interaction energy (F^^) seems to 
be a robust indicator of covalent interaction energy. Atomic promotion energies are attributed 
to the bonds (E^'^'") to derive Hartree-Fock bond energies (E^^*^). 



B. Consistent decomposition of molecular properties 



TablefV] displays the kinetic energies from the single-index atomic density matrices (see Eq. 



( 16 1) calculated with the matrix elements of the finite Aug-cc-pVTZ basis set. They are compared 



to the atom-condensed kinetic energies t^ [Eq. (36)] and t^' [Eq. (37)], obtained by numerical 
computation in r-space of the molecular kinetic energy integral weighted by the Hirshfeld-I atomic 
weight functions. We confirmed that all entries in table |V] are converged with respect to grid size. 
Note that for the homonuclear diatomic case, all values coincide, since there is only one way to 
partition the molecular kinetic energy over equivalent atoms. 

the atom-condensed kinetic energies t'\ and t^' differ significantly 



As discussed in section 



IV 



(e.g. a difference of 0.3 Hartree for C in CO) pointing to the inherent ambiguity of combining 
atomic weight functions with the kinetic energy operator. These ambiguities are absent when the 
molecular IDM itself is partitioned. For infinite basis sets, values for the kinetic energy of the 
IDM fragments t^ should equal values for t^. It is clear from the table that for the present Aug- 
cc-pVTZ basis set significant differences still occur. For example, the difference \tQ — tc\ is more 
than 0.06 Hartree in CO. We checked for CO that this difference vanishes in the Aug-cc-pVQZ 
basis set i\t^-tc\= 0.002 Hartree). 
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TABLE V: Kinetic energy (t^) of the molecular Hirshfeld-I fragments compared to the Hirshfeld-I frag- 
mented molecular kinetic energies (t^) and (t^). Values are presented relative to the kinetic 
energy (t^) of the neutral isolated atom. Computations were performed at the ROHF/Aug-cc- 
pVTZ level of theory. 
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C. Considerations about computational efficiency 



In this section, the computational efficiency of the approach to calculate atomic self-energies 



and the interaction energies in one-electron Hilbert-space (see Eqs. (16) and (18) ) is compared 
with the computational cost of the approach to compute these quantities using numerical integra- 



tion in r-space (see Eqs. ( [36| ) and (|40|)). When relatively small molecules and basis set sizes 
are considered, the matrix approach is computationally much less demanding than the r-space 
approach. Figure [T] displays the ratio of the computational costs of both approaches as a function 
of the number of atoms. In order to obtain sufficiently accurate energy integrals (condensed to 
atoms and atom pairs) in the r-space approach, one must go to very large grids, making the cal- 
culations quite time consuming. On the other hand, in the matrix approach, moderate grids suffice 
to construct the atom condensed atomic overlap matrices and the kinetic energy integrals over the 
molecular orbital basis are easily computed analytically. 
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FIG. 1: Cafcufation times for the energy decomposition (at the ROHF/Aug-cc-pVTZ level of theory) using 
the matrix approach (ti) refative to the computational cost of using an r-space approach (^2) as a 
function of the number of atoms in the moiecuie. 



Going to larger molecules and very large basis sets (e.g. Aug-cc-pVQZ), at one point the 
matrix approach will get computationally more demanding (at least at the Hartree Fock level) than 
the approach of Eq. ( [36] ). In the former method all virtual molecular orbitals are included in the 
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calculation, while in the latter method one only needs the occupied molecular orbitals. However, 
the crossover point is far from being reached in the present set of molecules. 

VIII. CONCLUSIONS 

In this paper, we applied a method for decomposing the one-electron molecular density matrix 
over the atoms in the molecule to calculate atom-condensed energy contributions at the Hartree- 
Fock level. 

In our opinion, this method is preferable to other approaches because it determines (Hartree- 
Fock) energy terms naturally associated with molecular fragments, whereas most other approaches 
fragment the molecular energy without assuring that there is an underlying electronic structure 
from which these energy fragments can be calculated explicitly. Since the Hartree-Fock energy 
can be expressed directly in terms of the molecular IDM, energy terms are derived from the 
molecular IDM fragments. For the current study, the IDM fragments are consistent with the local 
and positive definite Hirshfeld-I partitioning of the electron density. We have shown that without 
this mathematical consistency, the results are necessarily ambiguous. 

Unlike in most cases where numerical integration is used intensively, the new scheme requires 
only the Hamiltonian matrix elements (one- and two-electron integrals) in Hilbert space, that are 
routinely computed in ab initio programs, and the atomic overlap integrals. Analysis of the com- 
putational efficiency of the new approach versus the more common one based on numerical inte- 
gration of the molecular energy integrals weighted by atomic weight functions, shows that there 
is a trade-off between on the one hand the number of integrals that need to computed and the 
size of the grids required to reach an acceptable level of accuracy. When using numerical integra- 
tion for, e.g., the kinetic energy, relatively fewer integrals need to computed but these are found 
to require large integration grids. When using the new density matrix formulation, more atomic 
overlap integrals are required but these can be computed with sufficient accuracy already for mod- 
erate grids. Analysis of the computational efficiency shows that the density matrix approach is 
computationally much more efficient. 
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X. APPENDIX 



When the molecular density matrix p(r,r') = J2i d'iiJi{r)4'i{r') is partitioned, the atom- 
condensed kinetic energy based on pAir, r') can be written as: 

1 



tA 



drdr'5{r — r') V'^pA{r, r' 



~2 J drdr'5{r — r') 
~2 J drdr'5{r — r') 



1 



V p(r, r')-[wA{r) + WA{r')] 



(VV(r,r')) ■ ]^[wA{r) + WA{r')] 



+p(r,r') ■ -V^WA{r) + (Vp(r,r')) • VwA{r) 



J2di dr \tPi{r) {V^tPi{r)) WA{r) 



1 



+ (Mr)? 2 i^^^Air)) + Mr) (VV^,(r)) ■ VwA{r) 



(46) 



Removing V'^WA{r) and VwA{r) using partial integration, Eqs. (33) and (36) are shown to be 
equivalent: 

tA = -l^di I drwA{r) ilJi{r)V^tlJi{r) 



1 
2 

1 
2 



~9 J drwA{r) [Mr) (V'^.(r)) 

drdr'5{r — r')wA{r) W'^pir, r') 



(47) 



To demonstrate that the atom-condensed kinetic energy based on pA{r,r') does not depend 
on the representation of the kinetic energy operator, we repeat the former exercise starting from a 
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different respresentation of the kinetic energy operator: 



tA 



drdr'6{r - r') V ■ VpA{r, r') 



^ j drdr'6{r — r') 
- J drdr'5{r — r') 



V ■Vp{r,r')-[wA{r)+WA{r')] 



(V ■ VV(r, r')) ■ -[wA{r) + WA{r')] 



+ (Vp(r, r')) ■ ]^V'wA{r') + (V'p(r, r')) ■ ]^VwA{r) 
^^rf, /"rfr \{ViJi{r)fwA{r) + iJi{r){ViJi{r))-VwA{r) 



(48) 



Removing Vwa(?') using partial integration, Eqs. (34) and (36) are seen to be equivalent as well: 

tA = ^^rf^ j drwA{r) [{Vij,{r)f - V ■ {UrW^{r))] 



drdr'5{r — r')wA{r) \/'^p{r, r') 



.h 
Ia- 
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